﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using DotSpatial.Topology.KDTree;

namespace HAMath
{
    public class OridinateKriging : Kriging
    {

        public OridinateKriging(List<KrigingPoint> Pts)
            : base(Pts)
        {
            //this.Points = Pts;
            //this.Tree = this.ConstructKdTree(Pts);
            this.KrigType = KrigingType.Ordinary;

        }


        /// <summary>
        /// 查找符合参数距离范围内的点
        /// </summary>
        /// <param name="SearchDistance"></param>
        /// <returns></returns>
        private List<KrigingPoint> GetNearstPt(KrigingPoint pt)
        {
            List<KrigingPoint> NearstPts = new List<KrigingPoint>();

            //第一步：先搜索符合搜索选项最大个数的离散点
            double[] xy = new double[2];
            xy[0] = pt.X;
            xy[1] = pt.Y;

            object[] Nearst = this.Tree.Nearest(xy, Math.Min(this.SearchOption.MaximumNum, this.Points.Count));  //有可能离散点的个数比选项设置时的最大个数要小，所以当查找选项设置的最大个数是出错

            foreach (object obj in Nearst)
            {
                NearstPts.Add(obj as KrigingPoint);
            }

            //第二步：根据搜索半径查找半径范围内的离散点
            for (int i = 0; i < NearstPts.Count; i++)
            {
                double distance = MathHleper.GetDistance(pt, Points[i]);
                if (distance > this.Ellipsoid.SearchRadius) NearstPts.Remove(NearstPts[i]);
            }

            //第三步：再判断点的个数是否符合搜索选项最小个数
            //这一步要用到KdTree结构                       

            if (NearstPts.Count < this.SearchOption.MinimumNum)
            {
                DotSpatial.Topology.KDTree.KdTree temTree = this.ConstructKdTree(NearstPts);
                object[] NearstObj = temTree.Nearest(xy, this.SearchOption.MinimumNum);
                NearstPts = new List<KrigingPoint>();
                foreach (object obj in NearstObj)
                {
                    NearstPts.Add(obj as KrigingPoint);
                }

            }

            return NearstPts;
        }


        /// <summary>
        /// 对某一点进行插值，得到插值
        /// </summary>
        /// <param name="pt"></param>
        /// <returns></returns>
        public double GetValue(KrigingPoint pt)
        {
            //第一步：求待求点的邻近点
            List<KrigingPoint> nearstPts = this.GetNearstPt(pt);

            if (nearstPts.Count == 0) return -999;
            
            double[] weight = new double[nearstPts.Count + 1];
           
            if (this.Model.Range > 0 && this.Model .C >0)  //变程不能为0。拱高也不能为0，为了少报错，这里也按变程为0时的方法走
            {
                //第二步：求待求点与邻近点模型变差值
                double[] matrix1 = new double[nearstPts.Count + 1];

                for (int i = 0; i < nearstPts.Count; i++)
                {
                    double h = MathHleper.GetDistance(pt, nearstPts[i]);
                    matrix1[i] = this.Model.CalModelVariation(h);
                }
                matrix1[nearstPts.Count] = 1;


                //第三步：求邻近点各点之间的模型变差值
                double[,] matrix2 = new double[nearstPts.Count + 1, nearstPts.Count + 1];

                for (int i = 0; i < nearstPts.Count; i++)
                {
                    for (int j = i; j < nearstPts.Count; j++)
                    {
                        double h = MathHleper.GetDistance(nearstPts[i], nearstPts[j]);
                        matrix2[i, j] = this.Model.CalModelVariation(h);
                        matrix2[j, i] = matrix2[i, j];
                    }
                }

                for (int i = 0; i < nearstPts.Count + 1; i++)
                {
                    matrix2[i, nearstPts.Count] = 1;
                    matrix2[nearstPts.Count, i] = 1;
                }

                matrix2[nearstPts.Count, nearstPts.Count] = 0;

                //第四步：解克里金方程，求各点权重。

                double[,] InverseMatrix = MathHleper.GetInverseMatrix(matrix2);
                weight = MathHleper.MultiplyMatrix(InverseMatrix, matrix1);
            }
            else
            {
                //当变程为0时，各个点的权重一样，所以其实就是用平均法求平均值
                for (int i = 0; i < weight.Length; i++)
                {
                    weight[i] = 1.0d / weight.Length;
                }
            }
            //下面这句的主法是通过矩阵相除的方法来求权重，速度快些
            //weight = MathHleper.DivideMatrix(matrix1, matrix2);

            //第五步：求待求点的值
            double result = 0;

            for (int i = 0; i < nearstPts.Count; i++)
            {
                result = result + weight[i] * (double)nearstPts[i].Value;
            }

            return result;

        }




    }
}
